%% ========================================================================
%          Emissions trading system
%            
%          Main file for structural analysis.  
%
%%=========================================================================
% Change to location of master matlab
% Set up
run('setup.m');

%% ========================================================================
%                       Data: Read plant emissions and bids
%%=========================================================================

% Data files
filename_plants = [ data 'structural_model_input.txt'];
filename_bids   = [ data 'bids.txt'];
filename_period = [ data 'period_statistics.txt' ];

% Set rng. Note that rng also set at top of etsCounter.run method,
% ets.etsAbatementModelsWrapper, and ets.estimateControlStringency so that
% things can be run in different orders without results changing from
% running straight through.
rng(0);

% Read in plant and bid data
edat = etsData( filename_plants, filename_bids, filename_period, false);

display(edat.period);

% Initialize estimation object
ets = etsEstim( edat.plant, edat.bid );
% Estimate All ETS abatement models in treatment (in ets.fits)
ets = ets.etsAbatementModelsWrapper();
% Estimate command & control stringency in control
ets = ets.estimateControlStringency;


%%% Include if you want to only estimate results using pre-covid data
edat_precovid = etsData( filename_plants, filename_bids, filename_period, true);
ets_precovid = etsEstim( edat_precovid.plant, edat_precovid.bid );
% Estimate All ETS abatement models in treatment (in ets.fits)
ets_precovid = ets_precovid.etsAbatementModelsWrapper();
% Estimate command & control stringency in control
ets_precovid = ets_precovid.estimateControlStringency;

%% ========================================================================
%                       Counterfactuals: with real data estimates
%%=========================================================================
% Here we run 4 different counterfactuals: 
% 1) Isoelastic MC, 2) Isoelastic MC with heterogeneity by equipment,
% 3) step MC no free cyclone, 4) step MC with free cyclone

period_num = 4;
cap = [ 2e4 : 1e4 : 3e5 ]; % Calculate on the basis of caps imposed

% Initialize counterfactual object with estimated model
ctr = etsCounter( ets, cap );

%% Isoelastic cost function
ctr.scenario.costs = 'iso';

% Without plant heterogeneity
ctr.scenario.iso_type = 'no_hetero';
ctr.scenario.regime = 'Market';
ctrMIso = ctr.run;  
ctr.scenario.regime = 'Command';
ctr.scenario.ermodel = 'heat_eps';
ctrCIso = ctr.run;
% writout table for trading panels
temp = unique(ctrCIso.plants(:,["id_plant","id_gpcb"]));
temp = renamevars(temp, "id_gpcb", "gpcb_id");
writetable(outerjoin(temp, ctrCIso.output_itk,  ...
    Type="right", ...
    MergeKeys=true),...
    strcat(data, "/cc_heateps_output_itk_gpcb_id.csv"));
% Evaluate abatement cost elasticity at control emissions
ctrCIso.abateCostElasticityControlEmissions( period_num )

% With plant heterogeneity
ctr.scenario.iso_type = 'hetero';
ctr.scenario.regime = 'Market';
ctrMIsoHet = ctr.run;  
ctr.scenario.regime = 'Command';
ctr.scenario.ermodel = 'heat_eps';
ctrCIsoHet = ctr.run;

%% Step cost function
ctr.scenario.costs = 'step';
% set no heterogeneity for clarity, but doesn't impact step function
ctr.scenario.iso_type = 'no_hetero';

% Without free cyclone
ctr.free_cyc_step = false;
ctr.scenario.regime = 'Market';
ctrMStep = ctr.run;
ctr.scenario.regime = 'Command';
ctr.scenario.ermodel = 'heat_eps';
ctrCStep = ctr.run;


% With free cyclone
ctr.free_cyc_step = true;
ctr.scenario.regime = 'Market';
ctrMStepFree = ctr.run;  
ctr.scenario.regime = 'Command';
ctr.scenario.ermodel = 'heat_eps';
ctrCStepFree = ctr.run;

%% Plots using ctr Objects
% iso elastic cost figures
figures_path = figures_iso_fullnohet;
ctrPlotsWrapper(ctrMIso, ctrCIso, ets, edat, figures_path, optfig)

% iso elastic with heterogeneity cost figures
figures_path = figures_iso_fullhet;
ctrPlotsWrapper(ctrMIsoHet, ctrCIsoHet, ets, edat, figures_path, optfig)

% Step function no free cyclone figures
figures_path = figures_stepnofree;
ctrPlotsWrapper(ctrMStep, ctrCStep, ets, edat, figures_path, optfig, true, ctrMIso)

% Step function free cyclone figures
figures_path = figures_stepfree;
ctrPlotsWrapper(ctrMStepFree, ctrCStepFree, ets, edat, figures_path, optfig, true, ctrMIso)



%% Table of emissions and cost outcomes across all periods
table_path = tables;

% Iso-elastic cost function
ctrAbatementCostTableWrapper( ets, edat, period_num, table_path, '/Table_D1_A.tex', 'iso', 'no_hetero', false)
% Precovid iso
ctrAbatementCostTableWrapper( ets_precovid, edat_precovid, period_num, table_path, '/Table_F7_A.tex', 'iso', 'no_hetero', false)
% Iso-elastic cost function with heterogeneity
ctrAbatementCostTableWrapper( ets, edat, period_num, table_path, '/Table_E1.tex', 'iso', 'hetero', false)
% Precovid iso with heterogeneity
ctrAbatementCostTableWrapper( ets_precovid, edat_precovid, period_num, table_path, '/NULL.tex', 'iso', 'hetero', false)
% Step function with free cyclone cost function
ctrAbatementCostTableWrapper( ets, edat, period_num, table_path, '/Table_D1_B.tex', 'step', 'no_hetero', true)
% Pre Covid Step function with free cyclone cost function
ctrAbatementCostTableWrapper( ets_precovid, edat_precovid, period_num, table_path, '/Table_F7_B.tex', 'step', 'no_hetero', true)

%% Make benefit cost analysis table
run("cost_benefit_calculation/benefit_analysis_calculation.m")
